### to replicate Table 2 of Cao and Ward: Transnational Climate Governance Networks and Domestic Regulatory Action; 
### for a special issue in International Interactions;
### Thursday, 2.June.2016
### XC at Penn State: email xuc11@psu.edu for questions. 

dd <- c("C:/XunCao/Research/TCG/final/data code/table 2")  ### you need to set up your own working directory though. 
setwd(dd)

set.seed(300)

options(scipen=200)
library(xtable)

### 1: run the model using gbme.r from Peter Hoff: 
Yd<-dget("dat.tcg2")[[1]] ## this is the matrix for the dependent variable which counts the number of shared TCGs between countries;
Y<-Yd
n<-dim(Yd)[1]             ## you need to specify n to run the gbme.r
Xd<-dget("dat.tcg2")[[2]] ## this is N by N by 15 an array for dyadic covariates; 


source("gbme.r")          ## source the code for gbme.r;
gbme(Y=Y,Xd=Xd,fam="poisson",k=3,odens=100, directed=F, NS=200000, oround=6)       # run; this will generate two output files ==> this is going to take a while ... maybe more than one hour;
                                                                                   # you need OUT for the posterior distributions; 


    
      
### 2: summarize the results: after running the MCMCs;
out<-read.table("OUT",header=T)             #read in output
burn<-500 # the length of burn in: 500 maybe is a little too much ... but we have 2000 rows in OUT anyways; 
tit<-c("trade proxy", "distance", "domestic legislation", "energy production",
                                             "green NGOs",
                                             "trade share",
                                             "real gdp pc", 
         "polity", 
         "cci", 
#         "oilgascoal", 
         "CO2 pc", 
         "EU", 
         "bank deposits per gdp", 
         "federation", "population", "iso14000 per gdp")#, "country random")

## plot figure 6:
par(mfrow=c(3,5))
for (i in c(1:15)){plot(density(out[-c(1:burn),i+3]), main=tit[i], cex.main=1.7); abline(v=0, lwd=2); abline(v=quantile(out[-c(1:burn),i+3], c(.025, .975)), lty=2)}


## we need some summary statistics of the posteriors:
Sum.post<-NULL
for (j in c(19, 4:18, 20:dim(out)[2], 3)){
                               sum.post<-c(quantile(out[,j][-c(1:burn)], c(.025, .975))[1], mean(out[,j][-c(1:burn)]), 
                               quantile(out[,j][-c(1:burn)], c(.025, .975))[2] , 
                               sd(out[,j][-c(1:burn)])                               )
                               Sum.post<-rbind(Sum.post, sum.post)
                               }

Sum.post<-round(Sum.post, digits=4)
colnames(Sum.post)<-c("2.5%", "Mean", "97.5%")
rownames(Sum.post)<-c("constant", tit, colnames(out)[c(20:dim(out)[2], 3)])
Sum.post
xtable(Sum.post, digits=c(rep(4,4)))
